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ABSTRACT 


Validated results are presented for the new 3D body of revolution finite element- 
boundary integral code. As usual, a Fourier Series expansion of the vector electric and 
magnetic fields is employed to reduce the dimensionality of the system and the exact 
boundary condition is employed to terminate the finite element mesh. The mesh 
termination boundary is chosen such that it leads to convolutional boundary operators for 
low O(n) memory demand. Improvements of this code are discussed along with the 
proposed formulation for a full 3D implementation of the finite element-boundary integral 
method in conjunction with a CGFFT solution. 

OBJECTIVE 

The objective of this task is to develop innovative techniques and related software 
for scattering by three dimensional composite structures. The proposed analysis is a hybrid 
finite element-boundary integral method formulated to have an O(n) memory demand. This 
low storage is achieved by employing the FFT to evaluate all boundary integrals and by 
resorting to an iterative solution algorithm. Particular emphasis in this task is the 
generation of software applicable to airborne vehicles and the validation of these by 
comparison with measured and other reference data. Because the approach is new, a step 
by step development procedure has been proposed over a three-year period. During the 
first year the technique was developed and implemented for two-dimensional composite 
structures. Support software for the two-dimensional analysis such as pre- and post- 
processor routines were developed during the second year and a formulation was also 
developed and implemented for three-dimensional bodies of revolution. Finally, during the 
third year, we will develop, implement, and test the method for arbitrary three dimensional 
structures. 
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BACKGROUND 

Interest in three-dimensional (3-D) methods has increased in recent years, however, 
the associated demands in computation time and storage are often prohibitive for electrically 
large 3-D bodies. Vector and concurrent (i.e. hypercube, connection, etc.) computers are 
beginning to alleviate the first of these demands, but a minimization of the storage 
requirements is essential for treating large structures. 

The traditional Conjugate Gradient Fast Fourier Transform (CGFFT) method [ 1 ] - 
[4] is one such frequency domain solution approach which requires O(n) storage for the 
solution on n equations. This method involves the use of FFTs whose dimension equals 
that of the structure under consideration [5] - [7] and, therefore, demands excessive 
computation time when used in an iterative algorithm. Also, the standard CGFFT requires 
uniform rectangular gridding that unnecessarily includes the impenetrable portions of the 
scatterer. With these issues in mind, a new solution approach is proposed for solving 
scattering problems. The proposed method will be referred to as the Finite Element- 
Conjugate Gradient Fast Fourier Transform (FE-CGFFT) method. 

During last year's effort the FE-CGFFT method was developed for two- 
dimensional scatterers where the finite element mesh was terminated at a rectangular box. 
Inside the box boundaries, Helmholtz equation is solved via the finite element method and 
the boundary constraint is obtained by an appropriate integral equation which implicitly 
satisfies the radiation condition. Along the parallel sides of the box, this integral becomes a 
convolution and is, therefore, amenable to evaluation via the FFT. The dimension of the 
required FFT in this hybrid method is one less than the dimensionality of the structure thus, 
making it attractive for 3-D simulations. Also, because it incorporates the finite element 
method, the FE-CGFFT formulation remains valid regardless of the structure's geometry 
and material composition. 

The proposed method described in the University of Michigan Report 02592 1-6-T 
(see also [8]) is similar to the moment method version developed by Jin [9]. Jin's method 
was in turn based on work published in the early 70's by McDonald and Wexler [10] who 
introduced an approach to solve unbounded field problems. The proposed method is also 
similar to other methods (a few of which will be mentioned here), neither of which 
provides a storage reduction comparable to the proposed FE-CGFFT method. The 
unimoment method [11] uses finite elements inside a fictitious circular boundary and an 
eigenfunction expansion to represent the field in the external region. The coefficients of the 
expansion are then determined by enforcing field continuity at the finite element (FE) mesh 
boundary. The coupled finite element-boundary element method [12] uses the finite 
element method within the boundary and the boundary element method to provide the 
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additional constraint at the termination of the mesh. Unlike the proposed method, the 
solution in [12] was accomplished by direct matrix inversion (as in [9]), and the outer mesh 
boundary is not rectangular to take advantage of the FFT for the evaluation of the boundary 
integrals. 


EROGRESS 

Part of our efforts in this task were devoted to debugging and validating the three 
dimensional body of revolution (BOR) code developed in the previous months. The 
analysis associated with this code is described in the U of M technical report 025921-1 8-T 
where we also include validation data obtained over the past two months. Some of these 
are shown in figures 1-3 and refer to an ogive, a circular cylinder and a sphere. 
Unfortunately, it was found that as the bodies became larger the system's condition 
deteriorates and this was traced to the pulse basis formulation employed for the 
discretization of the boundary. Through several tests we have now shown that A 
Galerkin’s linear basis formulation will correct the convergence difficulties. For example, 
this formulation was already employed in solving large systems (with more than 120,000 
unknowns) associated with the scattering by frequency selective surfaces (FSS) and large 
Mates. As shown in figure 4-5, the Galerkin's formulation with linear basis permitted a 
solution of this system in less than 70 iterations! In comparison, the pulse basis-point 
matching formulation would require several thousand iterations before reaching 
convergence. Consequently, we are in the process of incorporating the Galerkin's linear 
basis formulation into our existing 2D and 3D BOR codes. Further, it was found that the 
ech area converges much sooner than the mean square error and permitted us to speed-up 
solution time. 

During this last quarter we also began the development of the proposed finite 
element formulation for general non-symmetric inhomogeneous bodies. The basic discrete 
elements in this case are tetrahedra in conjunction with edge-based expansion functions. 
The associated finite element formulations is described in Appendix A and we are now in 
the process of implementing it. Initially, the finite element mesh will be terminated by a 
fictitious absorbing layer whose dielectric parameters were determined by a minimization of 
the reflection coefficient over the entire range of incidence angles. For a three layer 
coating, each of thickness 0.05 wavelengths, it was found that their respective dielectric 
properties to minimize the reflection coefficient over all angles of incidence are 
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e r i=(-0. 1249205, -1.731605), p rl =(- 1.03 1792; -0.1039932) 
e r 2=(0.040699530, 0.1750280), ^n=(0.3155941, 0.3190330) 
e r 3 =(- 1.278644, 0.9625375), ^ r3 =(-0. 17213 15, -5.389832) 

The corresponding plot of the reflection coefficient as a function of incidence is given in 
figure 7 along with scattering patterns based on the proposed termination model. As seen, 
for the chosen fictitious absorber the reflection coefficient is less than one percent for 9 up 
to 62 degrees and less than 2 percent for 0 up to 77 degrees. For the same error criteria, 
the corresponding angles associated with the second order Pade ABC are 35 and 41 
degrees, respectively. The fictitious ABC has, therefore, a substantially better performance 
over the existing ABCs, and its effectiveness will be examined further in the next few 
months. 

The three dimensional finite element meshes required in the analysis will be 
generated by SDRC IDEAS and we have already began to develop the software for 
transforming the output of this commercial package to the input files of our analysis codes. 
Similar drivers were already developed for the two dimensional code which was developed 
last year. 

Finally, during this period we performed extensive testing of the two-dimensional 
code and have in the progress developed several new pre- and post- processing algorithms 
for this code. Two of the new geometries (see fig. 8 and 1 1) whose scattering was 
computed with our 2D finite element - CGFFT code are displayed in figures 9, 1 1 and 12. 
These represent airfoil configurations, one of which is coated with a dielectric material. 

CONCLUSIONS 

The project continues to evolve in accordance with our original plan and schedule. 
Most importantly, so far, our expectation of the finite element CGFFT formulation have 
been realized and we are, therefore, pleased with its performance for the intended 
applications. 

TRANSITIONS 

All of our efforts in the next six months will be devoted to the development of the 
3D finite element boundary integral code for arbitrary structures. In the immediate future 
we will also pursue improvements for our existing codes primarily directed at speeding the 
convergence of the CG or BiCG algorithm. 
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Body of Revolution enclosed in a fictitious finite length cylinder. Area between 
the Scatterer and the fictitious cylinder is discretized via the finite element method. 
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Figure 1: TM and TE bistatic scattering 
pattern from a perfectly conducting sphere of 
radius for axial incidence. 
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Figure 2: TM and TE bistatic scattering pattern from a perfectly conducting circular 
cylinder of length 1A. and radius 0.1 A. for axial incidence, (a) modes 0-2, (b) converged 
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Figure 3: TM and TE bistatic scattering pattern from a perfectly 
conducting ogive length IX and maximum radius of 0.088X for axial incidence. 
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Figure 4: Scattering by a 12 x 12 FSS array; comparison of the exact solution 
(solid line) with an approximate result obtained by truncating the infinite FSS. 
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Figure 5: Scattering by a 24 x 24 array; comparison of the exact solution 
(solid line) with an approximate result obtained by truncating the infinite FSS. 
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Figure 7: Evaluation of the new Fictitious ABC (a) geometry (b) reflection coefficient 

plot (c) H z Field on a PEC Circular Cylinder (d) E z Field at a distance X/10 from the 
surface of a PEC Circular Cylinder. 
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Figure 8: Geometry and finite element mesh of the illustrated coated trailing edge 
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Figure 9: E and H polarization scattering patterns for the configuration shown 
8. Comparison of results between the FE-CGFFT and the moment method. 
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Perfectly conducting airfoil 

All dimensions in wavelengths. The airfoil section is made by 5 arcs: 



• OA : straight line 

• Afi : circle of radios R2 = 7Ao and of center 02 

• BC : polynomial parametric equation 

• CD : polynomial parametric equation 

• DO : circle of radius Al = 9 Aq and of center 01 
The polynomial equation are given by: 

*(«) = *(«) * EjLifc • a 1 *"* 1 0 < « < 1 

^ and k for BC arc a,- and 4t for CD are 

«i * 4.61149 kt = 1.53278 o t * -1.62131 4i = -0.12563 

«t s -12.11403 4, = -3.22680 e, * 4.54389 6,= 0.30612 

o» * 8.88606 4,* 0.90615 a, = -3.32901 4, = 0.16113 

04 * 0.44275 4< * 1.24623 o« = 0.35663 4 4 = -0.00343 

a»* -0.51440 4,* 0.40927 a« = -1.40984 4,* -0.51216 

a*= -6.86720 4,* -0.41007 a* = -5.40756 4» = -0.23609 

Figure 10: Geometry of a PEC Airfoil whose scattering is given in figures 11 & 12. 
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Appendix A 


Finite element formulation for tetrahedral 
elements and edge-based expansion basis 

1 Derivation of finite element equations 

Let us consider a three dimensional inhomogeneous body occupying the 
volume V. In order to discretize the electric field E inside the body, we 
subdivide the volume V into a number of small tetrahedra, each occupying 
volume V € (e = 1, 2, • • • , M ) with M being the total number of tetrahedral 
elements. Within each tetrahedron, the electric field satisfies the vector 
wave equation 

V x —V x E — A^e r E — 0 (1) 

Mr 

where fi T is the permeability of the medium, e r is the medium permittivity 
and k 0 is the free space wave number. The next step is to expand the 
electric field within V e as 

E = ££'W‘ (2) 

J = 1 

where W* are edge-based vector basis functions and E * denote the 
expansion coefficients of the basis, all defined within the volume V e . W* is 
tangential to the jth edge of the eth tetrahedron with zero tangential 
component along the other edges of the tetrahedral element. On 
substituting (2) into (1), we obtain 

£i;fvxivx w ;-^w')=o (3) 

j=\ \ rr ) 

In order to solve for the unknown expansion coefficients Ej, we take the dot 
product of (3) with W' and then integrate the resulting equation over the 
element volume V t (Galerkin’s technique). The wave equation thus reduces 
to 
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J2 Ej f v w? • (v x J -' V x W' - ^rW’j du = 0 (4) 

The first term in the integral of the above expression can be simplified by 
using basic vector identities. Since 

W'-fvx-VxW;) = V- — (VxW')xW, e +-(VxW‘).(VxW]) 

V t*r J f*T J f^r 

the divergence theorem can be readily applied to (4) resulting in the 
following expression: 

0 = £ E ) { X (“( V X W * ) • (' V X W i) “ k ^r WJ.W^ dv 

+ £fc (VxW ' )xW -) ds } (5) 

where S e denotes the surface enclosing V e . Using vector identities , (5) can 
be further simplified to yield the weak form of Maxwell’s equation: 

[ (— (V x W t e ) . (V x W‘) - k 2 0 e r W'.W*) dv = ju>fi 0 <f W?.(n x H )ds 

j = l \ftr / ''Se 


where n x H is the tangential magnetic field on the exterior dielectric 
surface. Equation (6) can be conveniently written in matrix form as 

m[E*] = m 


a-, = [ f— (v X W‘) . (V X w;i - ) dv (S) 

JV*. \/^r / 

B\ = jup 0 jf w , e .( n x H)ds (9) 

On assembling all the M tetrahedral elements that make up the geometry, 
we obtain a system of equations whose solution yields the field components 
over the entire body. Therefore, summing over all M elements, we have 
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( 10 ) 


which gives 


M M 

= £[B1 

6—1 6=1 

[A][E] = [B] 


( 11 ) 


where [A] is a N x N matrix with N being the total number of edges 
resulting from the subdivision of the body and \E\ is a N x 1 column vector 
denoting the edge fields. Due to the continuity of the tangential component 
of the magnetic field at the interface between two dielectrics, an element 
face lying inside the body does not contribute to [5] since the surface 
integrals over the faces of adjacent tetrahedra cancel each other. As a result, 
\B ] is a column vector containing the tangential magnetic field only over 
the exterior surface of the body. Equation ( 11) can therefore be written as 

A SS E, + A„E, = H* 

A, S E S + A„E, = 0 (12) 


where the subscript s denotes the edges on the surface and i represents the 
edges inside the body. It is thus readily seen that (11) relates the electric 
field inside and on the surface of the body to the on-surface tangential 
magnetic field. 


2 Basis functions 


Vector fields within tetrahedral domains in three dimensional space can be 
conveniently represented by expansion functions that are linear in the 
spatial variables and have either zero divergence or zero curl. The basis 
functions defined below are associated with the six edges of the tetrahedron 
and have zero divergence and constant curl. Assuming the four nodes and 
the six edges of a tetrahedron are numbered according to Table 1 , the 
vector basis functions associated with the (7 — i)th edge of the tetrahedron 
are defined as 


W 


7— t 


fj—i + g 7 _, x r, r in the tetrahedron 
0, otherwise 


(13) 


where i = 1 , 2, • • • , 6 and f and g are constant vectors. On direct 
evaluation, it is readily seen that 
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V-W, = 0 

VxW, = 2g, 


(14) 

(15) 


Since the complex scalar Ej in ( 2) is the projection of the electric field 
onto the jth edge of the tetrahedral element, 

(16) 


on jth edge ~ 


where is the Kronecker delta. Solving ( 13) and ( 16) for the unknown 
vectors yield [1] 


fV-i 


g7 


h-i 
6V ril 
6^67— 

6V 


x r, 


(17) 

(18) 


where V is the volume of the tetrahedral element, e, = (r,- 2 — r tl )/6 t - is the 
unit vector of the ith edge and 6, = |rj 2 — r„ | is the length of the ith edge. 
All distances are measured with respect to the origin. 

Since there are two numbering systems, local and global, a unique global 
direction is defined (e.g., always pointing from the smaller node number to 
the larger node number) to ensure the continuity of n x E across all edges. 
This implies that ( 13) should be multiplied by (-1) if the local edge vector 
(as defined in Table 1) does not have the same direction as the global edge 
direction. Even though Wj forces no conditions on the normal component 
of E, it has been shown[2] that the continuity of electric flux can be satisfied 
within the degree of approximation with the above formulation. Finally, 
since V • W, = 0 the electric field obtained through ( 2) exactly satisfies the 
divergence equation within the element, i.e. V • E = 0. Therefore, the finite 
element solution is free from contamination of spurious solutions [2], 


3 Mesh termination 

Differential equation methods, such as finite elements, can only solve 
boundary value problems. Since electromagnetic problems are open 
boundary-infinite domain types, a means to truncate the solution domain 
to lie within a finite boundary must be found. On this boundary, a 
condition is enforced thus ensuring that the fields will obey the Sommerfeld 
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radiation condition at distances asymptotically far from the object. These 
absorbing boundary conditions (ABCs) have a significant advantage over 
the global methods of solving unbounded problems using finite elements in 
that they are local in nature. Due to this, the sparse matrix structure of 
the finite element formulation is retained. One disadvantage, however, is 
that ABCs are approximate and do not model the exterior field exactly. 

The objective of absorbing boundary conditions is to truncate the finite 
element mesh with boundary conditions that cause minimum reflections of 
an outgoing wave. These ABCs should provide small, acceptable errors 
while minimising the distance from the object of interest to the outer 
boundary. This minimal distance is required to reduce the number of 
unknowns in the problem for computational efficiency. A three dimensional 
vector boundary condition will be investigated here for terminating the 
finite element mesh of the body described in section 1.1. We begin with the 
Wilcox representation^] of the electric field which has an expansion 


E(r) 


e -ikr 


r 


f. A„(M) 

' ip 71 

n= 0 ' 


(19) 


From ( 19), we get 

V x E = jjifcf x 1 ^ lE 


” nA n< 

r 2 * r n 

7 71 = 1 ' 


( 20 ) 


where A nt = r x A n is the transverse component of A n and, for a vector F, 
DjF is given by 


DiF = 


sin6 




1 

'dF r 

nnjn/] h 1 ^ 

6 + 


1 

sin6 

89 

— siTiv r 

F m 


( 21 ) 


Using the recursion relation 

-2jknA nt = n(n - l)A n _ lt t + Z^An-! 


where 



(D At + D, A.) e + (d a* + d, a.) i 

8.4^ 1 , 2cosOBA* 

dO sin 2 6 n sin 2 0 80 

J_dK _ _1_ t 2coseaAj 

sinO d<j> sin 2 0 n sin 2 0 d<f> 

and D is Beltrami’s operator [3], we can derive the representation correct to 
t~ 4 . Applying the recursion relation in ( 20) yields the desired relationship 
for the vector ABC: 


D 4 A n = 
Dg A n = 
D$ A n = 


V x E = a(r)E + 0{r)D 4 E 


( 22 ) 


where 


c*(r) 

^(r) 



1 1 


2 jkr 2 (1 + 1/jkr) 


(23) 

(24) 


The ABC formulated above is applicable to spherical boundaries and hence 
would be storage intensive and numerically inefficient when used to 
terminate the mesh of long and thin geometries. It would be highly 
desirable to choose an outer boundary that conforms to the shape of the 
object. An approximate boundary condition based on the asymptotic 
representation of fields for a two dimensional scalar problem has already 
been derived [4]. It is the author’s intention to extend the derivation of the 
two dimensional scalar boundary condition to a three dimensional vector 
absorbing boundary condition for an arbitrary outer boundary. 


4 Solution of the finite element equations 

An inspection of ( 11) reveals that for an inhomogeneous body, there is no 
a priori information about the tangential magnetic field over the exterior 
surface of the body. Relation (11) therefore contains two unknown vectors, 
[E\ and [B], and thus another condition is required involving the two 
variables to permit an evaluation of the fields inside and on the surface of 
the body. This condition relating the tangential electric field to the 
tangential magnetic field on the surface is provided by ( 22). Since the 
ABC in ( 22) refers to the scattered field, we can rewrite it as 
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V x E 


: - a(r)E; + /?(r)D.E: 

H; = -i- [a(r)E; + 0(r)B,E;] 

i Jfi 

= fCEl (25) 


where K = ^ [a(r) + /?(r)D 4 ] and the subscript s denotes the field on the 
surface and the superscript s represents the scattered field. Since the total 
field is a sum of the incident field and the scattered field, therefore from 
( 25), we obtain 


h; = ice* 

H, — H' s nc = K (E 3 — E‘ nc ) 


(26) 


Substituting ( 26) into (12) and simplifying gives 

- K)E a + A ai Ei = H* nc -/CE* nc 

.A,’ s E a + j4 t iEi = 0 


(27) 


The above equation can thus be solved for the unknown electric fields both 
inside and on the surface of the body. 
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TABLE I 

Tetrahedron Edge Definition 
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